An experimental test of the local fluctuation theorem 
in chains of weakly interacting Anosov systems 
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Abstract 

An experimental test of a large fluctuation theorem is performed on a chain 
of coupled "cat maps" . Our interest is focused on the behavior of a subsystem 
of this chain. A local entropy creation rate is defined and we show that the local 



fluctuation theorem derived in |G1] is experimentally observable already for small 
subsystems, i.e. the finite size effects, for which the theory says little, are small 
even in small subsystems of the systems considered. 

1 Introduction 

A general law governing the fluctuations of the phase space contraction rate has been 
derived under various chaoticity assumptions, | GC1 ], | GC2 |, [MR|, Q, | LS1 ], | LS2 |, 
iQ, a review is in |G3], and it has been observed in several experiments, [BGG], 

|LLg, |cg, |Bcg, ilrbi . 

Consider a chaotic {i.e. Anosov) dynamical system S : x ^ Sx, on a phase 
space Ai of "size V" which we shall call "volume" and which will be the number of 
interacting subsystems composing our "spatially extended" system, see Sec. 2 for a 
precise definition. Let /i be the natural stationary state, i.e. the SRB distribution, 
and suppose that 5 is a time-reversible dynamics. One can study the fluctuations of 
the phase space contraction rate, per unit volume, averaged over a time r: 
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where r/_|_ is the infinite time average of ry(a;): = limT-^oo "J")^^ X]J=o ^('5-'^)- 
The value r/_|_ will be called the "average phase space contraction^^ per unit time and 
volume. It will also be called the average entropy creation rate. 

The probability distribution Tfrip) of the fluctuations p in the stationary state 
verifies, |]G1|, for large r: 
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T Vr]_^_ p + 0{V) "fluctuation theorem". 
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The size V is denoted 0(1) in all the papers to which we refer because the authors 
keep the number of degrees of freedom fixed but here it will be a parameter. 



The normalization proportional to V, used in (1.1), takes into account that in the 
models that we shall consider dissipation will occur in a translation invariant way so 
that T] will have a value typically proportional to V, hence an average Vf]_^ with r]_^ 
essentially independent ofV, i.e. of the spatial dimensions of the system (for V large). 
This implies that we expect, both from a physical and from a mathematical point 
of view, that in macroscopic systems the observable p has practically unobservable 
fluctuations on the scale involved in (|l.2| ). 

However if one looks at the motion of a subsystem localized in a spatially micro- 
scopic region V and if one considers a local entropy creation rate then fluctuations 
are far more probable. Hence observable consequences of a local fluctuation theorem 
are likely to become possible. 

It is not clear which would be, in any particular case, the appropriate definition 
of local entropy creation rate, see |BCL]. The simplest situation seems to arise in the 



case of a chain of V (V large) very chaotic systems weakly coupled to each other by a 
translation invariant coupling (see below). For such model systems there is a rather 
natural definition of local entropy creation, or local phase space contraction: at least 
it became natural to us after the first few attempts at performing the numerical study 
that we present in this paper on the basis of other definitions, that in retrospect were 
not as natural and which failed to give results of any interest. 

The definition of local entropy creation rate for the systems that we consider here 



has been proposed and discussed in [Gl|, where a local version of the fluctuation law 
has been derived. _ 

Consider a system with a "large" phase space Ai = 'S>J=o -^o where Mq is the 
phase space of a "single" system in a periodic chain of length V: hence a point in Ai 
representing a configuration of our chain will be a IZ-ple of coordinates (xq, • • • , Xy_-^). 
Let F be a small fixed interval [0,V — 1] C [0,1^ — 1], and let r]{x) be the "local 
phase space contraction rate" (yet to be defined and in general different from the 
r/(x) above) in the subsystem described by the coordinates (xq, . . . , xv~i). Define 
r/+ = limT-_+oo SJ=o vi'S-'x), and extend the definition ( p!^ ) above to y < y as 



E viS^^). (1.3) 



Then, if tTt{p) denotes the probability distribution of the above p, it verifies, see [Gl|: 



log =Tyr?+j9 + o(T + y). (1.4) 

where one should note that t + V \s (half) the size of the boundary of the space- 



time region involved in the definition oi p in (1.3). When V = V the error is only 
0{V) rather than 0(r + V), see (|1.2[ ), because there is no "vertical" boundary and, 
therefore, no boundary effect related to it, see |G1]. In ( |1.4| ) and in the cases that we 



shall consider, r]^ can be replaced by rj^ because it will he r]+ = r]_^ + 0{V ^) 
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The relations (|1.2| ), ( |1.4D are the objects of our experimental test in a concrete, 
simple, model for A^o and S. 

One can ask "why to make the test if there is a theorem that predicts the results?" 
This is true: but the theorem only leads to estimates of the remainders in ( |1.2D , ( |1.4| ). 
Such estimates are not quantitatively derived in the proof, because various constants 
in the bounds are only shown to be finite with no attempt at computing them. 
They can be, however, computed following a rather general approach (see |G1[| where 
the general theory of [PS], | BK1 |, [ BK2| ], |JP] is used in an essential way) and, as 
usual when applying general mathematical theorems to concrete cases, we expect the 
estimates to be very poor. 

Therefore it makes sense performing an experiment to check if the predictions 
of the theorem are actually observable: a situation similar to the one met in the 
theory of phase transitions and in the numerical experiments performed on the two 
dimensional Ising model. The latter model is exactly soluble in the sense that many 
of its properties can be analytically computed: nevertheless experiments on its a 
priori known properties make sense and are important to test the feasibility and the 
reliability of simulations on systems with an enormous number of degrees of freedom 
and to test the accuracy and reliability of methods that one really wants to apply to 
more complicated systems not amenable to an exact solution. 

The global fluctuation laws ( |l.2D are universal in many systems where, however, 
they are extremely difficult to observe: in this paper we give evidence that their 
local counterparts can be observed quite easily because they manifest themselves 
observably already in rather small subsystems. 

This is quite analogous to the theory of density fluctuation in a gas: one cannot 
observe density fluctuations on large gas samples because their likelihood decreases 
exponentially with the volume. Nevertheless such fluctuations are of great interest 
(their value will give us the compressibility of the gas) and they can be measured by 
observing them in small samples and by making use of the fact that the logarithms 
of the distribution of their averages over a given volume V are proportional to it so 
that once the fluctuations are known in a small sample one can infer their (otherwise 
practically unobservable) probability in a large sample. 

In section 2 we describe the particular system considered in our experiment and 
we define the local phase space contraction rate following [Gl|; in section 3 we expose 
the experimental results of our simulations. A few comments will follow in section 4. 
An appendix discusses briefly our method to estimate the experimental errors that 
are characteristic of our simulations. 



2 Description of the system 

Let M' be a two-dimensional torus, M' = [0, 1] x [0, 1], and x a pair of coordinates 
on it. A Arnold cat map S' : M' M' is the dynamics on M' generating from an 
initial datum x = a trajectory t ^ x* 

x^+^ = Sx^=(} 2)^* modi. (2.1) 
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Consider a finite number V of copies of such systems and label each copy with an 
index i G {0, . . . ,V — 1} identifying Xq and xy (periodic boundary conditions). 

Let A4' = Mj- be the phase space of such new "extended" system, let x = 

'S>iQV—i ^ point in it and let S' = <S>i(zv '^i ("unperturbed") evolution: this 
means that each coordinate Xi evolves under the map S' into SiXi = (S'x)i. We now 
consider a perturbation of the map S' 

S{x) = S'{x) + eJ^{x) (2.2) 

where is a "nearest neighbor coupling" so defined that: 

S{x)- = S-{xi) + e/.(xi_i,Xi,Xi+i). (2.3) 

We choose /. = (/«, 0) with /j very simple. Namely, if Xj = {Cj,S,j)j for j = 0, . . . , F — 
1, and if we set aj = 27r(^j_i + + ^j+i) and Pj = 27r(Cj-i + Cj + Cj+i): 

fi = ^[sm{ai + Pi) + sin p,] (2.4) 

On the boundaries we simply set = Cy_i; = ^v-i^ ~ ''O' ~ ^O' 
periodic boundary conditions. 

The system {M',S) is, for small perturbations (i.e. small e but for all V), an 
Anosov map by the remarkable extension of structural stability to chains of Anosov 
maps with short range interaction developed in |PS|. However the theory of the 
fluctuations that we test here applies directly only to ^'reversible systems, see [pCl 



GC2| and below. 



A map 5 is a time reversible map for a time reversal operation T if X is a smooth 
isometric involution I : M < — > A4, l'^ = Identity, such that IS = S^^I. 

Therefore the system that we shall consider in the experiment cannot just be the 
one described so far because the map S is not time reversible (at least not for any 
simple map I). We shall, however, consider a closely related system which has a 
time-reversal symmetry and which we now describe. 

Let M = M' M' and define: 

S{x,y) = {S[x),S-\y)), (2.5) 

The system (A4,5) is (trivially) reversible because the map I{x,y) = {y,x) is a time 
reversal for it (see []G2| ). 

Denote x = {x,y) a point in the phase space Ai: the ''phase space contraction 
rate" of {Ai,S) is defined, see [ GC1 |, in terms of the determinant Ay (x) = | det dS{x)\ 
of the jacobian dS{x) of the map S and of the determinant Ay{x) = \ detdS{x)\ of 
the map S. 

The global entropy creation rate for this system is defined as minus the logarithm 
of the Jacobian determinant of S, A{x) = — | det dS{x)\, |GC2f| , rj{x) = — log A(x): 

fjy(x) = — log Ay (x) = — log I det 55(2;) I 

/dS{x) 
= -log det /^~.~ 1 (2-6) 

= -log(Ay(x)(Ay(5y))"'' 
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If fij- are the forward and backward statistics of the volume measure {i.e. the SRB 
distributions for S and S~^) the perturbation makes, for generic perturbations /, see 



| BGM ], the system {Ai,S) dissipative, in the sense that: 

- / ^i±{dx) log A±H^) = V±y > 0. (2.7) 

at least for e ^ small. This is a property that we checked in our particular case 
(the form in |2.4| was so chosen to make the algebraic part of this check simpler) and 
which is a quite general fact, see [Q. 

Note that time reversal symmetry implies that r}_^_ = r]_ and r]^ = 2rj ii rjV is the 
average of — log (x) . 

Fixed a subvolume V, i.e. a subinterval of V, let dSv{x) be a block matrix 
of the Jacobian dS{x) corresponding to the coordinates Xy of this subsystem. Set 



Av{x) = \detdSvix)\ and define the local phase space contraction rate, see |G1[| , as 



riv{x) = -log (^Avix)[AviS-^y)) (2. 



that will also be called the local entropy creation rate in the volume V (we imagine 
V <^ V). We will study the fluctuations of the averages over a time r of r]y(x) and 
of riv{x), see (0), (0). 

3 Experimental results 

We considered chains of V = 8, 16, 32, 64 maps as defined in Sec. 2, and for each V 
subsystems of size V = 1,3,6 consisting of V consecutive elements of the chain. The 
perturbation parameter e was given the values 0.05, 0.10, 0.15 but the complete tests 
were done only for the value 0.10 and the other values were used only for the purpose 
of checking/testing the codes. 

According to the definition in Sec. 2 we associate with each chain site two pairs 
of coordinates each defining a point on the 2-torus: the collection of the first V 
pairs evolves under the action of a forward map S, while the collection of the second 
pairs evolves under S~^. The inversion of S has to be done numerically, because 
although the map S is very simple, its inverse cannot be elementarily expressed: we 



used a general Newton's algorithm (called Newton-Raphson algorithm in |NR]) which 
strongly limits the maximal size of the perturbation e but which, for the relatively 
small values of V we considered, is competitive with the 0(y)-size algorithm that we 
designed specifically for our problem. 

The total number of degrees of freedom is 4:V, hence it went up to 256 although 
the results that we report deal with the V = 16, i.e. with 64 degrees of freedom cases. 

For such values of e the Lyapunov exponents can be checked to be still very close 
to the ones corresponding to the case e = 0. Just as an example we give in Figure |l| 
the graph of the 64 Lyapunov exponents corresponding to a model with V = 16 and 
e = 0.10. 

Note that the intrinsic complexity of the inversion of the map S is not the only 
obstacle to the inversion of the map: as e increases the system becomes more and 
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Figure 1: Plot of the Lyapunov exponents for the system with V = 16, e = 0.10. The 
64 exponents (each "cat" has 2 components that evolve forward and 2 that evolve 
backwards in time) are ordered in couples A+ , A_ coupling the greatest positive with 
the greatest negative, the second greatest positive and the second greatest negative 
and so on. 



more dissipative and at some critical value stops being conjugated to the unperturbed 
system. For instance the attractor ceases to be dense in phase space, and density is 
a condition for the validity of the fluctuation theorem, see |GC2|. Although this new 
situation is interesting we shall not study it here, see [BG|. 

Our simulations were carried out over 4 • 10^ -4- 300 • 10^ iterations^, depending on 
the system dimension and on the perturbation parameter e0. 

The figures |2[ ^, |5| refer to the analysis of data produced by the system with 

V = 16, e = 0.10 for = 3 • 10*^ iterations. The values of F = 8,32,64 as weh 
as e = 0.05, 0.15 have also been considered and we do not report the results. We 
only say that they were very similar, with the obvious variation that in the cases 

V = 32, 64 substantially less statistics was available for the analysis of the global 
fluctuations (which become too rare as the size of the system increases). 

Our data and error analysis follows from the one described and applied in [ BGG ] , 
see the appendix below for further details. 

Once the values of e, V and V <V are fixed we measure the fluctuation of entropy 
creation rate p (see ( |1 . 1\) or ( |1.3| )) averaged over the iteration time r ranging from 2 
to 400, discarding tq iterations between each measure to let them decorrelate. We set 

^To let each simulation take few days of CPU time on a modern PC. 

■^The number of steps the algorithm takes to calculate the backward evolution depends on e: i.e. 
for e = 0.05 the algorithm converges (in the sense that the precision reaches the numerical double 
precision accuracy) in 5 or 6 steps while for e = 0.15 we have to use a separate algorithm to search 
for a "good" starting point for the Newton algorithm to converge in 50 steps. 
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Figure 2: Histogram of the distribution 7r(p) for the chain of = 16 subsystems and 
sub-chains of F = 1,3,6 with e = 0.10 and r = 40. The simulation was carried 
over N = 10* iterations. The amplitude of the distribution increases in smaller 
subsystems: large fluctuations are more probable and the fluctuation theorem can be 
easily observed. The probability of observing negative values of p in the full system 
is so small that we could obtain just few experimental points. The error bars (barely 
visible in the scale of the figure)are discussed in the appendix. 



To = 10 because we know from the computation of the Lyapunov exponents, reported 
in Figure that the maximal and the minimal ones are such that tq = 10 is about ten 
times the characteristic time of the system^. For the same reason we discard the first 
10 iterations in order to reduce the effects of transients while the motion approaches 
the attractor. 

Then we compute the empirical probability distribution of fluctuations of the 
phase space contraction rate Tfriv) for the whole chain and TTrip) for the subsystem. 
Figure 

■^The exponents are of the same order of the unperturbed positive exponents that are exactly 
log(i(3 + ^)) -0.98. 
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Figure 3: Plot of x{p), (|3J]), for the system with V = 16, e = 0.10, V = 1,3,6 and 
r = 40. The solid line is the linear least-square fit; the dashed line is the line with 
slope 1. In smaller subsystems large deviations are more probable (the tails of TTrip), 
Figure ^, go far below 0) and the fluctuation theorem prediction for r = +oo can be 
observed more easily than in the full system. 



To check the predictions of the global and local fluctuation theorems, ( |1.2[ ) and 
( |1.4| ), we verify that the quantities 

^(P)= ^l°g|T^' -(P) = ^>s^) (3-1) 

are linear in p, see Figure ^, with slope Xt- 

x{p) = XtP-, Xt ^ 1 for T ^ oo (3-2) 

We computed Xt with a linear least squares fit. Results, plotted in Figure ^ 
show that Xt does not contain 1 within the error bars only for small values of r. We 



attribute this to finite-size effects, | BCL |, i.e. as a manifestation of the 0(r~^) and, 
respectively, 0(r~^ + V~^) of the corrections in ( |1.2D , ( [1.41 ). We also checked that 
rjj^ = rjj^ + 0{y~^) and in fact the correction in 0(y~^) is not detectable for V > 16. 
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Figure 4: Plot of the slope, Xt, of x{p) for the system with V = 16, e = 0.10 and 
Vo = 1, 3, 6. In smaller subsystems the fluctuation theorem holds for larger ranges of 
values of r. For r small, Xt does not contain 1 within the error bars: we attribute 
this to finite-size (of r) effects. The straight line has the theoretical slope 1. 



A least square fit of the probability distributions W-r (p) and (p) gives a Gaussian 
as best fit (in the sense of | BGG |, [ GG |). 

If p was Gaussian then it would also verify (|3.2| ) with some slope (not necessarily 
1): however there is no reason for p to he a Gaussian variable, and in fact it can only 
be such for p rather close to its average value p = 1. In fact around p = 1 it looks 
Gaussian because of a central limit theorem generally valid for the SRB distribution 
of Anosov systems, see [Q, or of weakly interacting chains of such systems, see |PS|, 

To test whether p can be really taken Gaussian, we computed the second moments 
of the distributions vr^(p) and tTt{p)- Assuming that Tfrip) and TTrip) are Gaussian and 
agree with the central limit theorem (which, however, only really holds for |p — 1| < 
0((ry)~2) implies that their 2"*^ moments m2(r) and m2(T) should be linear in 
{tV)^"^ (to leading order in tV) and we would denote them m2(r) = [tV)^"^ A, 
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Figure 5: Plots of ^ for the system with V = 16, e = 0.10 and Vq = 1,3,6. The 
experimental values of AA/A differ from 0. The experimental value of A differs 
from the value Aq that would follow from a Gaussian assumption combined with the 
fluctuation theorem by ~ 1%, quite clearly visible beyond the error bars. The range 
of r in the graphs is over the values of r corresponding to Xr = 1) see Figure 



m2(T) = (tV) ^ A. The distributions of p would then be proportional to: 

e-rF{p-i)V(2A)^ ^-rV{p-if/{2A) ^3^3^ 

for some A, A respectively. 

Therefore the fluctuation relations ( |1.2| ) , ( |1.4D or ( |3.2| ) will imply that there should 
be a simple relation between the slopes 1 of x{p) and x{p) and the mean square 
deviations coefficients rn2(r) = {tV)~^A, m2(T) = {tV)~'^A of the random variable 
{p — 1), namely: 

2 _ 2 

— = A, — = A. (3.4) 

V+ V+ 

Calling Aq,Aq the values in ( ^ ) we deduce from the experimental of m2(r) 
and m2(r) as functions of Vt or of Vt, the actual experimental values A and A 
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(performing a least square fit). We then compare these values with the ones predicted 
in {i.e. with Aq = 2r?+\ = 2r?+^). 

AA A-Ao AA A-Ao 

-A- = TT' = — ^'-'^ 

In Figure ^ we plot AA/A only, of course, for the values of t for which the fluctuation 
theorem can he checked {i.e. for the values of r for which there is enough statistics 
to have significant results in spite of the error bars, see Figure §). The results of this 
test is plotted in Figure ^: The amplitude of the error bars doesn't let us see that the 
fluctuations of the entropy production are not gaussian (in a numerical test sense) 
although the value of is not zero. The error analysis seems to play a foundamental 
role in this test. 



4 Conclusions 

(1) The "fluctuation theorem", for our models, is an exaci result only asymptotically 
when r ^ oo in the global case and when r, 1/ — > oo in the local case. It appears, 
nevertheless, to be observable already for small values of r, V provided the "correct" 
definition of local entropy creation rate is applied (see (p.7|), finding it was one of the 
starting points of the present work, see |G1| ). The observability of the predictions of 
the theorem on small systems is evident for y = 16 as can be seen in Figure |2[ ^ 
^ We have performed similar experiments for systems with 1/ = 8, 32, 64 obtaining 
quantitatively similar results with (of course) lower statistics in the last two cases: we 
do not report them for brevity. We also tested, although not in the same detail and 
the same accuracy, the values of e = 0.05 and e = 0.15, the latter being substantially 
more time consuming: the results seemed to agree with the ones that we report. 

(2) Large fluctuations are practically unobservable in large systems: it's evident 
also in the cases we considered. The plot of the probability distributions vrT-(p) of 
the fluctuations of p in the whole chain has just few points for negative values of p 
(Figure ||) thus the fluctuation theorem can be checked for few values of r (Figure 

However the local fluctuation theorem holds for any small subsystem and can be 
checked for a larger number of values of r (Figure ^ even when the global theorem 
can not he verified. 

(3) The hypothesis of p being Gaussian for either the whole chain or for its sub- 
systems is theoretically "impossible" because large deviations of normalized sums 
of independent nongaussian variables {i.e. of averages of independent equally dis- 
tributed nongaussian random variables) are not Gaussian even in the most random 
systems like Bernoulli schemes, |FGP|| . 

If the distribution of p was a Gaussian centered at p = 1 the second moment of 
p — 1 should have the form iyT)~^A and the slope of the logarithm of the odd (in p) 
part of the distribution oi p should have the form 2Vt/A and it should he Aq = A, 
which can be interpreted as the validity of the Green-Kubo formula in the nonlinear 



regime we are considering, see |ES| and Sec. 5 in |BGG[| and |]G3| , [1G4 |. 
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We plan an experiment to investigate more accurately the results of Sec. 3 indi- 
cating that is significantly nonzero: this means that we are trying to find systems 
for which the numerical construction of the trajectories is less time consuming for the 
values of V that we consider. 

The nongaussian nature of the fluctuations has already been studied and clearly 
shown in the work [E£] at least in models out of equilibrium. Non gaussian behavior 



has been also studied and announced by |G[. We find it also in equilibrium (a result 
not reported here) as (of course) expected. 

(4) The above considerations are also supported by the fact that the nongaussian 
nature of the fluctuations is clearly demonstrated in a situation in which, neverthe- 
less, the linearity of the odd part of the logarithm of the distribution of p can be 
observed (although it is not clear whether it falls within the theory of reversible 
Anosov systems): this is a case quite different from the presently considered chains of 



Anosov systems, see |LLP], and the nongaussian nature of C(p) is there manifest for 
r relatively small because the distribution appears sensibly non symmetric around 
p = l. 



A Appendix: Error evaluation. 



We follow the error estimates described in [|GG |, | BGG |; we add here explicitly the 



description of the technique we use to estimate errors, which are not trivially due to 
statistical errors, on the distribution histogram, because in the present experiment it 
is the main error source. 

Let I{k) be the count of the number of values of p in the k-th interval [kd, kd + d] 
(with d > a fixed mesh) into which the p-axis is divided to perform the histogram, 
and let I'{k) = I(—k — 1). The main source of errors is that r/_|_ itself is affected by an 
error; the histogram error due to an uncertainty 6rj in rj^ is the error in the number 
count I{k) over the interval [kd, kd + d] given by 2^^/(/i;): because the error ^ = ^ 
in the measurement of p amounts to an uncertainty p^ on the apparent position of 

the two extremes of the interval [kd, kd + d] centered at p = {k + ^)d. 
So that, adding the statistical error 3y/T(k), the total error on I{k) is: 

p 5r, 



AI(k) = 3Jl(k) + 2 1(k) -^d (A.l) 
V d r]+ 

The corresponding variation of the observed value of x: 

7y+r\/ \k) 

has to be taken (recalling that the ideal result is x(p) = p = (k ^ \)(^'- 

1 , l{k) \ fbj] 3 3 ^ ^ \ /AON 

■log-F777T =P — + ^ + ^ + 2PT^ + 2p:7^ (A.3) 
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which we have added and/or propagated to the other errors (due to the large size of 
our statistical samples the two purely statistical contributions to dx turn out to be 
essentially negligible compared to the others). 
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